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ABSTRACT 

In this work we s tudy the properties o f the mass density field in the non-Gaussian world mod- 
els simulated by iGrossi et alj (l2007l) . In particular we focus on the one-point density prob- 
ability distribution function of the mass density field in non-Gausian models with quadratic 
non-linearities quantified by the usual parameter /nl- We find that the imprints of primordial 
non-Gaussianity are well preserved in the negative tail of the probability function during the 
evolution of the density perturbation. The effect is akeady noticeable at redshifts as large as 
4 and can be detected out to the present epoch. At z = we find that the fraction of the 
volume occupied by regions with underdensity S < —0.9, typical of voids, is about 1.3 per 
cent in the Gaussian case and increases to ^ 2.2 per cent if /nl = —1000 while decreases to 
~ 0.5 per cent if /nl = +1000. This result suggests that void-based statistics may provide 
a powerful method to detect non-Gaussianity even at low redshifts which is complementary 
to the measurements of the higher-order moments of the probability distribution function like 
the skewness or the kurtosis for which deviations from the Gaussian case are detected at the 
25-50 per cent level. 
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1 INTRODUCTION 

The renewed interest on studying deviations from Gaussianity in 
the primordial density fields has been recently prompted by new 
observations and theoretical progress. 

On the theoretical side, inflation-based non-G aussian (NG) 
mode ls have been thoroughly investigated (see, e.g.. lBartolo et al.l 
|2004 and references therein) and can be used to make robust pre- 
dictions for quantities that can be potentially observed lik e the mass 
functi on of virialized structures (Matarrese et al. 200Q), their bi- 
asing iMatar rese et al. 1986, ,: Grinstein & Wise 1986 ; Dalai et al., 
I2OO8 : Matarre se & Verdel l200a). the topology of the large scale 
structure (LSS) of the universe 1 Matsubarall2003l). and highe r-order 
clustering statistics like the bi-spectrum 1 Hikage et alj2006l) . 

On the observational side, the cosmic microwave background 
(CMB) has provided the most stringent constraints so far on the 
level of non-Gaussianity. The most r ecent analysis of the WMAP 
5-year temperature fluctuation maps jKomatsu et al.ll2008 ) shows 
that the primordial NG signal on the very large scales probed by 
CMB is tiny (see the discussion in SectO. Yet the latest claims 
of a positiv e detection of NG feature s in the WMAP 3-year data 
reported by lYadav & Wandel^ ( |2008|) has renewed the interest in 
NG models, shifting the focus to smaller scales where NG features 



can only be spotted through the analysis of the Large Scale Struc- 
ture (LSS) of the universe at much lower redshifts. In this respect, 
current observations and future planned experiments will certainly 
deliver new information capable of setting new and stringent con- 
straints on the amount of non-Gaussianity. 

In this framework numerical experiments represent the only 
tool capable of bridging the gap between theoretical analytic pre- 
dictions that usually involve simplifying hypotheses, and LSS- 
based observations, where non-linear effects play a non-negligible 
role. In this respect, different groups have performed new, in- 
dependent numerical experiments (see, e.g., Mathis et al. I I2OO4I : 
[kang et al. 2007; Grossi et al. 2007; Dalai et al. 2008). 

In this work we use the highest-resolution experiments per- 
formed so far of physi cally motivated NG models that have been 
described in detail by jOrossi et al.ll200l hereafter G07). In that 
work we have studied the evolution of massive dark matter haloes 
in non-Gaussian scenarios quantified by the mass function, and 
compared it with analytic predictions. In a subsequent paper we 
have used the same numerical experiments to investigate the topo- 
logical properties of t he mass density fi eld quantified by the 
Minkowski Functionals jHikage et al.ll2008l) . 

Here we focus on the probability distribution function of the 
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mass density field [PDF hereafter] which describes the probability 
that a randomly placed cell of specified volume and shape contains 
a specified density. In absence of primordial non-Gaussianity, the 
PDF has a Gaussian form initially, but develops a positive skew- 
ness and kurtosis at later times. Deviations from the original Gaus- 
sian shape originate from the non-linear evolution of density fluc- 
tuations and non-linear distortions induced by departures from the 
Hubble flow, if redshifts are used to estimate galaxy distances. In 
our analysis we will ignore the latter and use N-body simulations to 
account for non-Gaussianity induced by non-linear dynamics with 
the purpose of characterizing the best range of scales and density 
in which primordial non-Gaussianity dominates over late contribu- 
tions. 

Our approach should be regarded as complementary to the an- 
alytic methods that model the evolution of the PDF in a Gaussian 
scenario using either perturbation theory or excursion set methods 
(see, e.g.. Lam & Sheth 2008, and references therein). Assessing 
the goodness of these methods is beyond the scope of this work. 
However, as a by-product, we will compare our result to the log- 
normal model to evaluate whether this simple analytic prediction is 
accurate enough to allow measuring non-Gaussianity in the PDF or 
whether more sophisticated models are required. 

Finally, we point out that our analysis should be regarded as 
somewhat ideal since, like in majority of the theoretical PDF stud- 
ies, we focus on the effects on the non-linear evolution of density 
fluctuations ignoring not only redshift-space distortions but also de- 
viations from Gaussianity due to galaxy biasing. 

This paper is organized as follows. In Section|2]we introduce 
the NG model considered here. The characteristics of our N-body 
simulations are presented in Section |3] where we also give details 
on the method adopted to implement suitable NG initial conditions. 
In Section |4] we illustrate how primordial non-Gaussianitities af- 
fect the redshift evolution of different quantities like the cluster- 
ing properties (Section |4Tt , the probability distribution function of 
density fluctuations (Section|4j2} and its high-order moments (Sec- 
tion |4.3t We discuss the results and conclude in Section|5] 



2 NON-GAUSSIAN MODELS 

It is possible to model the level of primordial non-Gaussianity pre- 
dicted by a large class of models for the generation of the initial 
seeds for structure formation (including standard single-field and 
multi-field inflation, the curvaton and the inhomogeneous reheating 
scenarios) using a quadratic term in the Bardeen's gauge-invariant 
potential "I>, which, on scales much smaller than the Hubble radius, 
reduces to minus the usual peculiar gravitational potential: 

-1> = '1>l + /nl($l-{$l)) . (1) 

In the previous equation <I>l represents a Gaussian random field, 
while the specific value of the dimensionless non-lin earity param- 
eter /nl depends on the assumed scenario (see, e.g.. lBartolo et al.l 
[2004). Note that in the literature there are two conventions for equa- 
tion Q, b ased on the LSS and the CMB, respectively (see a dis- 
cussion m iMatarrese & Verdell2008l) . Here we follow the first one, 
for which $ is the gravitational potential linearly extrapolated to 
z = 0. 

It is worth stressing that detailed second-order calculations of 
the evolution of perturbations from the inflationary period to the 
present time show that equation l[T) is not generally valid. In fact 
the quadratic, NG contribution to the gravitational potential should 
be represented as a convolution with a kernel /nl (x, y) rather than 



a product (see, e.g.. iBartolo et al.ll2005l) . Since in this paper we 
will consider models for which |/nl| ^ 1, all space-dependent 
contributions to /nl can be safely neglected and /nl can be ap- 
proximated by a constant. In particular the simulations that we will 
analyse here assume /nl = ±100, ±500, ±1000. We notice that 
owing to the smallness of $, the contribution of non-Gaussianity 
implied by these values of /nl is always within the percent level 
of the total primordial gravitational potential, and does not appre- 
ciably affect the linear matter power spectrum. 

We notice that the considered range is larger than that al- 
lowe d by present CM B data. In fact, the analysis performed by 
Kom atsu et"ar] ( |2008|) on the WMAP 5-year maps finds 9 < /nl < 
111 for the 95 per cent confidence level. Our choice for the /nl 
range is motivated by two different reasons. First, the LSS pro- 
vides observational constraints which are a priori independent of 
the CMB. Second, /nl is not guaranteed to be scale independent, 
while the LSS and CMB probe different scales. Indeed some in- 
flationary sceiwi£s_do_2redi£HM scale-dependent /nl (see, 
e.g.. lChenll2005l : IShandera & Tveli2006l) . 



3 N-BODY SIMULATIONS 

The set of 7 simula tions presented here (see also lGrossi et alj2007l : 
iHikage et alj2008 h is characterized by different values of /nl, and 
adopts the same so-called 'concordance' ACDM model consisting 
of a flat universe dominated by a dark energy component given in 
the form of a cosmological constant A. The relevant cosmologi- 
cal parameters match those determined using the WMAP first-year 
data (Spergel et al. 2003), i.e. Q,mo ~ 0.3 for the matter density 
parameter, f^Ao ~ 0.7 for the A contribution to the density pa- 
rameter, h — 0.7 for the Hubble parameter (in units of 100 km 
s~^ Mpc~^). The initial power spectrum, which adopts a cold dark 
matter (CDM) transfer function, has a spectral index n — 1 and 
its amplitude guarantees that the r.m.s. of density fluctuations mea- 
sured in spheres of radius 8h^^ Mpc is as = 0.9. We recall that 
the recent analysis o f the WMAP 5-year data jKomatsu et al.l200^ ; 
iDunklev et alj|20()8l) provided a slightly different set of best-fitting 
cosmological parameters. In particular, they suggest smaller values 
for both (78 and QmO- Adopting the first set of parameters rather 
than the second does affect the outcome of some statistical analy- 
ses presented in this work. However, it does not influence the gen- 
eral results that focus on the relative differences between models 
characterized by a different degree of primordial non-Gaussianity. 

All N-body simulations have been performed at the CINECA 
Supercomputing Centre (Bologna) and at the 'Rechenzentrum der 
Max-Planck-Gesellschaft' (Garching) by using the publicly avail- 
able code GADGET-2 ( Springel 2005), which is a mo re flexible and 
efficient version of the numerical code GADGET dSpringel et al.l 
1200 ll) . thanks to the improved parallelization strategy. The gravi- 
tational acceleration is computed using the tree method, in combi- 
nation with a particle-mesh (PM) scheme for long-range force. In 
this work we follow the evolution of collisionless dark matter par- 
ticles, only. Therefore we are not making use of the hydrodynamic 
capability of the numerical code. 

In all experiments we loaded a computational box of (500 
Mpc//i)^ with 800"^ particles: the corresponding particle mass is 
then m — 2.033 x lQ^'^h~^ Mq. The gravitational force has a 
Plummer-equivalent softening length of e; = 12.5/i~^ kpc. The 
grid used by GADGET in the PM scheme was set to have 1024'' 
nodes. The runs produced 31 outputs from the initial redshift and 
the present time; the redshifts of the last 25 outputs have been fixed 



The density field in non-Gaussian models 3 



in order to have a comoving distance of 250/i~^ Mpc (i.e. half box) 
between them: this will allow to easily build past-light cone real- 
izations for different complementary studies already in progress. 

3.1 Setting non-Gaussian initial conditions 

The most critical step in performing the simulations is setting up 
suitable NG initial conditions without modifying the power spec- 
trum of density fluctuations, P{k), which has to be identical to that 
of the Gaussian case. 

To do that we start by performing a random realization of a 
Gaussian gravitational potential characterized by a power-law 
spectrum P{k) oc This is done in Fourier space using a 1024'' 
grid, by applying the standard technique based on the combination 
of a random phase with an amplitude extracted from a Rayleigh 
distribution corresponding to the desired P(k). This gravitational 
potential is then inverse-Fourier transformed back to real space to 
include the NG term: $nl = /nl ($l — As a final step 

we go back to fc-space to account for the CDM matter transfer func- 
tion. This procedure, wh ich improves the original one proposed by 
iMoscardini e"tai] ( ll99lh . guarantees that different NG models with 
different values of /nl have the same linear power spectra, all con- 
sistent with that of the Gaussian case, as we have checked. Besides 
the Gaussian case with /nl = 0, we have considered six different 
NG models characterized by /nl = ±100, ±500, ±1000. 

The initial conditions for our N-body simulations are then ob- 
tained by using the resulting gravitational potential to perturb a 
uniform distribution of particles according to the Zel'dovich ap- 
proximation. In order to avoid applying this technique beyond its 
range of validity, the initial redshift Zi is fixed in such a way that the 
maximum (absolute) value for the density fluctuations on the grid is 
smaller than 0.5. As a result the redshift at which we start the sim- 
ulation depends on the /nl value and varies in the range [73, 100] 
with the highest redshift corresponding to the largest /nl (abso- 
lute) value. In order to avoid spurious effects due to the non perfect 
isotropy of a grid, we prefer to perturb a glass-like distribution, ob- 
tained as a result of a N-body simul ation where th e sign of gravity 
has been reversed. As discussed by I White! jl994h , this technique 
avoids the presence of preferred directions in the cubic grid, so the 
initial distribution retains the uniformity feature. For this purpose, 
we replicate 8 times a 100'^ glass file along each Cartesian axis, 
exploiting the periodic conditions. 

In Fig.[T]we show the probability density function of the den- 
sity field, as measured in the earliest output (2 — 50) common to all 
our 7 simulations. It is computed on a 512^^ grid, using the Trian- 
gular Shape Cloud (TSC) scheme (see, e.g. lHocknev & Eastwood 
Il988l) . From the plot it is evident that the larger the value of /nl, 
the more pronounced is the high-density tail. This reflects on the 
values of the skewness of the density distribution {S'*), which 
ranges at 2: = 50 between 1.3 x 10"^ for /nl = —1000 and 
1.4 X 10""' for /nl = 1000. Note that the positive values of the 
skewness in models with negative /nl are produced by the first 
stages of the gravitational evolution. In fact, taking into account 
our sign choice in the relation between potential and density, we 
expect a negative (positive) skewness for the PDF when a nega- 
tive (positive) /nl is used. We have checked this by analyzing the 
simulated density field at the redshift corresponding to our initial 
conditions. We find that at Zi the skewness is —3.9 x 10~^ and 

1.2 X 10"^ for /nl = -1000 and /nl = 1000, respectively. 
The small asymmetry in the two values of the skewness is partially 
due to sample variance and also to the fact that the two quantities 
refer to different epoch since Zi = 96.7 for /nl = —1000 and 
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Figure 1. The PDF for the different models (as indicated in the plot) as 
computed at the earliest common redshift (2 = 50). 

Zi — 100.8 for /nl = ±1000. The redshift evolution of the proba- 
bility distribution function of density fluctuations and its moments 
will be discussed in detail in Sections [4.2l and l431 



4 RESULTS 

In this section we investigate the evolution of the dimensionless 
power spectrum A^(fc), and the one-point distribution function of 
the mass overdensity p{5, Rs), focusing on its high-order moments 
and low-density tails. 

4.1 Redshift evolution of clustering 

A first visual impression of how the value of /nl affects the clus- 
tering evolution can be gained by looking at Fig. [2] where we show 
maps of the mass density field at different redshifts ranging from 
2 = 5.16 and z — 0. Each slice has a side length of 500/i"^ Mpc 
and a thickness of 31.25/i~^ Mpc. The central column presents the 
results for the Gaussian simulation (i.e. /nl = 0), while the rela- 
tive residuals for /nl ~ —1000 and /nl = 1000 are shown in the 
corresponding panels of the left and right columns. The residuals 
are calculated at each pixel as 

Apa; = {px - po)/po , (2) 

where px represents the density for the NG map with /nl ~ 2; and 
po refers to the Gaussian field. All density fields in the plots have 
evolved from the same random realization of the underlying linear 
gravitational potential, irrespective of the statistics. For this reason, 
the maps of the residuals show similar features but with opposite 
signs. 

As a first quantitative statistical test, we consider the dimen- 
sionless power spectrum A^, defined as the contribution to the vari- 
ance of the fractional density per unit In k: 

- ^k - ^''''^'^ • 

In the previous formula P{k) is the matter density power spectrum. 

The effect of non-Gaussianity can be appreciated in Fig. |3] 
where we show A^ (fc) computed at four different redshifts, rang- 
ing from 2 = 3.96 to 2 = 0. For the sake of clarity we present the 
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Figure 2. Slice maps of the simulated mass density fields at 2 = 5.15, z = 3.96, z = 2.13 and 2 = (rows from top to bottom). The side length is 500/i^^ 
Mpc, while the thickness is 31.25/i~^ Mpc. The panels in the middle column show the logarithm of the projected density. The left and right panels are the 
relative residuals for the /NL=ilOOO runs (equation (2)). Each panel has the con'esponding colour bar; note that for the sake of clarity the considered ranges 
are different from redshift to redshift. 
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Figure 3. The redshift evolution of the dimensionless power spectrum A^. 
For clarity, only the non-linear scales are shown. Different lines refer to 
models with different primordial non-Gaussianity: /nl = (solid line), 
/nl = 1000 (dotted line), /nl = -1000 (dashed line). 



Figure 4. The logarithmic deviation R(k) = log A^j^^(A;) — log A^j^(fc) is 
shown for the same models and redshifts pres ented in Fig E] He re Ajj^(A;) 
is computed using the analytic expression of ISmith et al.l 
valid only for the Gaussian case. 



4.2 Probability distribution function of density fluctuations 



results only for the scales which are in the non-linear regime (i.e. 
for > 0.1 /i/Mpc). Moreover we show only the two most extreme 
NG models (/nl = ±1000), which are compared to the Gaussian 
reference one. 

As expected, at a given redshift, a positive value of /nl tends 
to increase the power at small scales, anticipating the formation 
of cosmic structures with respect to the Gaussian case. On the op- 
posite, negative /nl values delay structure formation. Fig.[3]shows 
that the effect is very small atz — O and increases with the redshift. 
However, even atz — 3.96 the largest deviations from the Gaussian 
case at fc = 0.5 h/Mpc do not exceed the 30 per cent level, i.e. well 
below those found in the negative tail of the PDF, as we will show in 
the next Section. In order to better quantify this effect, in Fig.|4]we 
present, for the same models and redshifts shown in Fig.|3] the loga- 
rithmic deviation, R(k) , of the meas ured A'^i^{k) from the analytic 
expression Ati-j(fc) that lSmith et al.l (2003) have obtained from the 
analysis of a large library of N-body simulations. We point out that 
this analytic model can be applied only to models with Gaussian 
initial conditions. Indeed in the simulation with /nl = we find 
R{k) = log A^i^(A;) - log At\(fc) < 0.01. In the two NG mod- 
els with /nl ~ ±1000 the dimensionl ess power spectrum deviates 
from the Gaussian case (and from the ISmith et^ ( |2003|) predic- 
tions) at fe ~ 1 /i/Mpc by ~ 0.04 at z = 3.96, \R\ ~ 0.05 
at z = 2.13 and z = 0.96, and \R\ ~ 0.03 at z = 0. Notice 
that at redshifts z > 4 (not shown in the plot), the differences be- 
tween the dimensionless power spectra computed for the Gaussian 
and NG models decrease because non-linear effects on the scale of 
interest become negligible at large redshifts. This result suggests 
that the best redshift range to spot primordial non-Gaussianity is 
2 < 2 < 4. 



The simplest way to characterize the cosinological density field is 
its one-point PDF, i.e. the probability of a density fluctuation 5 mea- 
sured at the generic position, on a suitable smoothing scale, Rs . If 
the PDF of the primordial density field is Gaussian then it remains 
Gaussian as long as density fluctuations evolve in the linear regime. 
When \5\ ~ 1, non-linear evolution strongly modifies the Gaus- 
sian shape of the original PDF. This picture has been confirmed by 
the results of different N-body s imulations perform ed in the cold 
dark matter framework (see, e.g. iKavo et al.ll200ll and references 
therein). These numerical experiments show that the evolved PDF 
is well approximated by a lognormal distribution: 



Pln{5) 



: exp 



[ln(l + S) + 
2S2 



1 + S 



(4) 



where the parameter S is related to the variance = {<5^) of 
the density fluctuation field: — lii(l ± a^). The validity of this 
result has been thoroughly checked in the weakly non-linear regime 
but also seems to apply well in the non-linear regime, at least for 
(J < 4 and S < 100. This result is further corroborated by the 
statistical analysis of the three-dimensional distribu tion of galaxies 
(see, e.g., ICroton et alj|2004l : iMarinoni et al.ll2005h . A discussion 
of the possible role of the lognormal distribution f or cosm ological 
structure formation can be found in Coles & Jonel ( Il99ll) . 

Here we extend the previous analysis by investigating to what 
extent the lognormal model can be successfully applied to describe 
the evolved PDF in scenarios with primordial non-Gaussianity. To 
do this we have computed the mass density field of our simula- 
tions using the TSC scheme to distribute each particle mass at the 
points of a regular grid. To study the clustering evolution we have 
considered different outputs corresponding to different redshifts 
(z = 3.96, 2.13, 0.96, 0). Moreover, since the non-linear scale in- 
creases with time we have also considered two different grids: a 



6 M. Grossi et al. 



finer one consisting of 512'^ points, corresponding to a smoothing 
radius of 7?s ~ 0.98 Mpc/h, suitable for the early stages of the evo- 
lution, and a coarser one, with 128^ points ~ 3.91 Mpc/h) for 
the low-redshift outputs. The PDF at the very early times {z — 50) 
has been already shown in Fig.[T] It is evident that deviations from 
Gaussianity are very tiny even for the most extreme NG models, as 
quantified by the low values of the skewness. 

In Fig. [5] we show the cumulative PDFs at four different 
epochs. We prefer to show the cumulative rather than the differen- 
tial PDF since, for a given overdensity threshold, the difference in 
the PDFs quantifies the difference in the volume occupied by struc- 
ture below threshold, i.e. a quantity that can be potentially mea- 
sured from real datasets. It is evident from the plots in Fig.[5]that 
the largest differences between models are in the low-density tail 
corresponding to log(l + S) < —1, as expected from a visual im- 
pression of Fig.[2]in which the most prominent features in the map 
of the residuals occur in regions that are below average density in 
the Gaussian case (shown in the central panels). These residuals 
are seen already al z = 5.16 and their magnitude increases with 
time. Negative /nl values increase the probability of underdense 
regions with respect to the Gaussian case. The opposite holds true 
for positive /nl. On the contrary, in the high-density tail, the differ- 
ences between the models are much less evident. The differences 
are more noticeable at 2 — 2.13 and decrease both at higher and 
lower redshifts, as already pointed out. 

Finally, departures from the Gaussian case are more evident 
when the density field is smoothed on a small scale (Rs ~ 0.98 
Mpc/h in the plots on the left hand side of Fig.|5]l, but can still be 
appreciated when smoothing on a scale of ~ 3.91 Mpc/h (shown 
in the panels on the right). Considering the density maps with 
Rs ~ 3.91 Mpc/h, commonly adopted in similar analyses, and 
an underdensity threshold S = —0.9, a valu e close to the mea n 
void underdensity in numerical simulations tColberg et alj|20()8h . 
we find that at z = only 1.3 per cent of the volume is below 
threshold in the Gaussian case. The difference with the NG models 
is striking, since the volume reduces by a factor 2.3 in the case of 
/nl = 1000 and increases by ^ 70 per cent when /nl ~ —1000. 
At higher redshifts the differences further increase. With the same 
overdensity threshold we find that if /nl = 1000 the volume oc- 
cupied by regions with S < 0.9 is a factor of ~ 6.5 smaller than 
in the Gaussian case while if /nl = —1000 the volume is ~ 17 
times larger. One should note, however, that the volume fraction 
below this density threshold rapidly decreases with the redshift. At 
z = 0.9 this fraction is only 0.007 per cent for /nl ~ 0. The 
fact that at all redshifts departures from Gaussianity preferentially 
occur in low-density environments strongly suggest that dedicated 
statistical tests related to the void probability could be successfully 
applied to spot primordial non-Gaussianity. It is worth pointing out 
that this is a promising observational test also at the present epoch 
since the volume occupied by voids, in which the NG features are 
still well preserved, increases with time. 

In order to check if a lognormal distribution provides a good 
description of the evolved PDFs even in the NG case, we com- 
pute the logarithmic deviation, AlogP = log PDF — logPLjv, 
where the reference lognormal Pln is computed from equation 
(|4j, adopting for the variance the value computed directly from 
the N-body outputs. We notice that the value of the mass variance 
that we measure in the Gaussian simulation agrees, at the few 
percen t leve l, with the theoretic al predictions of IPeacock & Doddsl 
"(Il996h andlSmithetalJ ( l2003h . hence confirming the results of 
lKavoetaljf200lh . 



The results are displayed in Fig. |6] for the same models, red- 
shifts and smoothing radii presented in Fig. [5] In general a log- 
normal distribution provides a reasonable fit to the evolved PDFs. 
The match is very good at intermediate overdensities (—0.2 < 
log(l + 5) < 0.6) and improves when the smoothing length in- 
creases. The deviation A log P increases towards both the high 
and the low-density tails. This trend is systematic and indepen- 
dent of the primordial non-Gaussianity. In particular the lognor- 
mal distribution tends to over(under)-estimate the probability of the 
under(over)-dense regions, and therefore cannot describe the statis- 
tical properties of both voids and highly-clustered regions. We re- 
call that analytic theoretical predictions for the PDF, different from 
the lognormal model, which agree at few pe r cent level in t he voi d 
regime have been very recently proposed bv lLam & ShethI ( l2008b . 

4.3 High-order moments 

A way to quantify the PDF properties is to measure its high-order 
moments. Rigorously speaking all moments have to be specified to 
fully describe a distribution function. However, for practical pur- 
poses a reasonable characterization of the bulk of the PDF can be 
achieved using the small-order moments up to the fourth. There- 
fore, here we consider the second, the third and the fourth-order 
moments of the distribution, namely its variance = (5^), skew- 
ness 7 = {S"^) and kurtosis k = {5"^). 

The gravity-driven clustering evolution amplifies these mo- 
ments. Their growth rate, however, does not depend on the level 
of primordial non-Gaussianity. This is evident from Fig. |7] where 
we show the redshift evolution of variance, skewness and kurto- 
sis relative to the PDF measured on the finer (512'^) grid, i.e. to 
a smoothing radius Rs ~ 0.98 Mpc/h. For the sake of clarity 
and to avoid overcrowding, we display the results for the Gaus- 
sian case (/nl = 0) and for the two most extreme NG models 
(,/nl ~ ±1000) only. The evolution of each moments is very sim- 
ilar in the three models, characterized by different line-styles. This 
is particularly evident for the variance, a'^, and reflect the fact that 
the power spectra of the various models remain almost identical 
throughout the evolution. Focusing on 7 we note that NG models 
with positive (negative) /nl, characterized by a positive (negative) 
primordial skewness, have values of 7 that are systematically larger 
(smaller) than the Gaussian case by about 25 per cent . This differ- 
ence, originated from the different initial conditions, is kept at the 
same level through the evolution history. Indeed, to a first approxi- 
mation, the gravitational evolution is expected to couple high-order 
moments to the variance (see below), the evolution of which has 
only a weak dependence on /nl, as witnessed by the fact that the 
evolution of P{k) depends only weakly on /nl (see Section l^Tt . 
This is also confirmed by the analysis of the other simulations of 
NG models with I/nlI < 1000 (not shown). These same consider- 
ations apply to the kurtosis, k for which, however, deviations from 
the Gaussian case at z = amount to ~ 50 per cent. 

Starting from Gaussian density fluctuations, perturbation the- 
ory shows that gravitational instability leads to hierarchical scal- 
ing relations between the variance of the PDF and its higher- 
order moments. For this reason it is useful to repeat the previ- 
ous analysis using two different quantities: the normalized skew- 
ness, defined as, 5* = 7/0-'* and the normalized kurtosis, K = 
(k — 3a*)/(T*'. Predictions of the perturbation theory by differ- 
ent authors hav e been successfully tested using N-body simula- 
tions (see, e.g.. iLucchin et al]|l994l and references therein). The 
results, however, revealed that the values of 5* and K depend on 
the shape of the primordial power spectrum and on the shape and 
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Figure 5. The redshift evolution of the cumulative PDF of the density field smoothed on two different scales. Panels on the left: smoothing radius Rs '^0.98 
Mpc/h. Panels on the right: smoothing radius Rs ~ 3.91 Mpc/h. The different panels refer to the four redshifts explored, indicated in the plots. Different 
line-styles refer to models with different primordial non-Gaussianity: /jml = (solid line), /nl = 1000 (dotted line), /nl = —1000 (dashed line). 




Figure 6. The logarithmic deviation of the PDF from a lognormal distribution, A log P, is shown for the same models and redshifts presented in Fig.|5] Results 
for smoothing radii Rs ^ 0.98 and Rs ~ 3.91 Mpc/h are displayed in the left and right panels, respectively. Different lines refer to models with different 
primordial non-Gaussianity: /nl = (solid line), /nl = 1000 (dotted line), /jml = —1000 (dashed line). 



scale R of the smoothing filter. For instance, with a Gaussian filter 
and a cold dark matter power s pectrum, typical values are 5* ~ 3 
and K ^ IG ( Kayo et alJlioOlh . These values also agree with the 
predict ions of the logno rmal model described by equation for 
which ( iKavo et alj|200lh 



S{R) =3 + a^{R) 
and 

K{R) = 16 + 15a^{R) + &o^{R) + g^{R) 



(5) 



(6) 



We recall that the hierarchical scaling relations are a predic- 
tion derived under the assumption of Gaussian fluctuations. For this 
reason, possible deviations from the theoretical predictions have to 
be regarded as a signature of primordial non-Gaus sianity. This is 
in esse nce the so-called skewness test suggested bv lColes & JonesI 
( Il99ll) and lSilk& Juszkiewiczl ( Il99lh . according to which depar- 
tures from a constant value S on scales in which fluctuations are 
still evolving in the weakly non-linear regime provide a power- 
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1 +z 

Figure 7. The redshift evolution of variance a'^ , skewness 7 and kurtosis 
K. Labels on the plot identifies the different moments. Different line-styles 
refer to models with different primordial non-Gaussianity: /nl = (solid 
line), /nl = 1000 (dotted line), /nl = -1000 (dashed line). 



ful test for NG models, as demonstrated by the analysis of first- 
generation NG simulations (Coles et al. 1993). 

We have repeated the same test using both the normalized 
skewness S and kurtosis K. The results are shown in Fig. [8] in 
which both quantities are shown as a function of the r.m.s. a. The 
different panels refer to different numerical simulations and com- 
pare the measured K (filled squares) and S (open circles) with the- 
oretical predictions of the lognormal model for which S{R) and 
K{R) are provided by equations ^ and ([6f (solid and dashed 
curves, respectively). Each panel combines results obtained at dif- 
ferent redshifts (z = 3.96, z = 2.13, z = 0.96 and 2 = 0) and us- 
ing different smoothing radii (Rs = (1, 2, 4, 8, 16, 32) x Ax, with 
At ~ 0.98 Mpc/h). The magnitude of both moments increases 
with a, independently of the value of /nl- Focusing on the top 
panel, we notice that in the low-cr regime (i.e. cr < 1) S{a) flat- 
tens but is not quite constant, which would be at variance with the 
predictions of the hierarchical model. However, we attribute this 
mismatch to a spurious effect deriving from the use of different 
smoothing scales. A similar trend is also seen in the two NG cases 
in the same Iow-ct range. We conclude that, given the smallness of 
the effect, it is difficult to disentangle the signature of primordial 
non-Gaussianity from the spurious signal produced by the varying 
smoothing radii, hence reducing the power of the skewness test. 
Using the normalized kurtosis has two potential advantages. First, 
the signal is more than one order of magnitude larger than the nor- 
malized skewness. And second, the differences between Gaussian 
and NG cases in the same a - range are significantly larger. How- 
ever, as shown bv lKavo etal] ( l200lh . kurtosis has larger error than 
skewness. Therefore to decide which statistics is more sensitive to 
primordial non-Gaussianity one should perform an accurate error 
analysis, which is beyond the scope of this work. Such analysis 
would also allow to estimate and correct for systematic effects in- 
duced by the finite size of the computational box although in our 
case, given the large size of the computational box, the correction 
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Figure 8. The normalized skewness S (open circles) and the normalized 
kurtosis K (filled squares) of the simulated density field as a function of 
u. The corresponding theoretical values predicted by equations jsj and (6) 
are shown with a solid and a dashed line, respectively. Different panels 
refer to models with different primordial non-Gaussianity. Top: Gaussian 
model /nl = 0. Middle: NG model with /nl = +1000 . Bottom: NG 
model with /nl = —1000 . Moments have been computed using differ- 
ent smoothing radii [Rs = (1, 2, 4, 8, 16, 32) X Ax, with Ax ~ 0.98 
Mpc/h] . The colours of the symbols correspond to measurements at differ- 
ent redshifts: z = 3.96 (red), 2 = 2.13 (blue), 2: = 0.96 (green) and z = 
(black). 

factor, estimated according to lKavo et al.l ( l200lh , is very small and 
can be safely neglected. 

The plot confirms that in the weakly non-linear regime a log- 
normal model predicts fairly well the normalized high-order mo- 
ments. Discrepancies become very significant during the non-linear 
evolution, i.e. when cr ^ 3 in the S case and when ct ^ 1 in the K 
case. We note that in the high-a regime the lognormal model sys- 
tematically over-predicts the amplitude of the normalized moments 
in both the Gaussian and the NG models. 



5 DISCUSSION AND CONCLUSIONS 

In this work we have investigated the properties of the mass den- 
sity field ex tracted from t he N- body experiments of NG scenarios 
presented in lGrossi et al. These models are physically mo- 

tivated in the context of the inflationary theory, and cover a large 
range of primordial non-Gaussianity: —1000 < /nl < 1000. 
Moreover the state-of-the art (mass, force and spatial) resolution 
of these simulations allows us to study the evolution of the pri- 
mordial non-Gaussianity on both linear and non-linear scales, from 
very high redshifts to the present epoch. 

Our results show that clustering properties are little affected 
by primordial non-Gaussianity, which is soon erased by the effect 
of non-linear evolution. Indeed, the power spectrum of our most 
extreme models matches the Gaussian case within 20-30 per cent. 
On the contrary, the one-point PDF of the mass density field rep- 
resents a promising tool to distinguish among competing NG mod- 
els. In fact, while in the high-density tail the primordial signal is 
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quickly obliterated by non-linearly induced non-Gaussianity, in the 
low-density regions the primordial non-Gaussianity stands out, ex- 
pecially at intermediate (1 ^ z < 2) redshifts. The differences are 
preserved up to the present epoch. 

This result suggests that statistics related to the void distri- 
bution could provide a powerful tool to detect primordial non- 
Gaussianity. Unfortunately, studying low-density environments is 
challenging from both a practical and a theoretical point of view. 
First of all, a theoretical difficulty is represented by the fact that no 
such statistical tool has been proposed in the context of NG models. 
Furthermore, even in Gaussian scenarios, there is no unique way of 
estimating a well-defined statistic like, for example, the void prob- 
ability function. In this respect, significant progress has been re- 
cently made to define efficient and reliable voi d finding algorithms 
(see for example the comparison presented in lColberg et al.ll2008l . 
and references therein). Finally, all statistics aimed at probing low- 
density regions are notoriously difficult to apply because of the 
sparseness of the data. However, this problem might be overcome 
by exploiting quasar high-resolution spectra thanks to the intrinsic 
relation between the properties of the void in the three-dimensional 
mass distribution and the "flux- voids" (defined as the connected re- 
gions in the one-dimensi onal flux distribution above the mean flux 
level; see IViel et al.l200g) . We plan to investigate this issue in detail 
by means of dedicated, high-resolution hydrodynamical simulation 
of NG models. 

Using spectra of distant objects to spot voids allows, in prin- 
ciple, to study their properties and therefore to investigate non- 
Gaussianity at all redshifts. At high z this method will comple- 
ment other statistics which are also effective in detecting NG 
features, like the occurrence of rare ev ents represented by the 
abundance of massive clusters (see, e.g.. iGrossi et aflboOTl) , the 
topology of t he ma ss distribution, that has been investigated by 
iHikage et alj j2008l) . Finally, in this paper we have studied the 
skewness and the kurtosis of the PDF and found that neither is 
very sensitive to primordial non-Gaussianity, even at very high 
redshifts. In this respect the bispectrum represents a much more 
powerful statistical tool that has been extensively used to con- 
strai n the amount of primordial non-Gaussian i ty in the CMB (see , 
e.g.,lKomatsu et al."2003'; 'Cabella et al."200^; 'Spergel et al.l l2007l : 
ICrem inelli et al. 2007; Komatsu et al. 2008). Also, it has been suc- 
cessfully applied to the existing galaxy surveys at lo w redshift to 
asses s the amount of non-linearity in galaxy biasing fVerde et al.] 
l2002l ; lGaztariaga & Scoccimarrdl2005l ; 1Pan & Szapudill2005h . In a 
future paper we plan to use our simulations to assess the power of 
the bispectrum in hi gh-redshift surveys as a pr obe of primordial 
non-Gaussianity (see lSefusatti & KomatsiJI 20071) . 
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